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Abstract 

A variety of system identification techniques are applied to the modelling of aircraft gas turbine dynamics. The motivation behind 
the study is to improve the efficiency and cost-effectiveness of system identification techniques currently used in the industry. Three 
system identification approaches are outlined in this paper. They are based upon: multisine testing and frequency-domain identifica¬ 
tion, time-varying models estimated using extended least squares with optimal smoothing, and multiobjective genetic programming to 
select model structure. © 2001 Elsevier Science Ltd. All rights reserved. 
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1. Introduction 

The modelling of gas turbines is a topic of great im¬ 
portance, given their widespread use in aero, marine and 
industrial applications. Three system identification tech¬ 
niques are employed in this paper, with specific applica¬ 
tion to an aero gas turbine engine. The techniques span 
the statistical and systems approaches described by 
Ljung (1996) in a recent survey of system identification. 
The motivation behind this work was to reduce test times 
while improving and quantifying model accuracy. The 
project involved three university teams, in collaboration 
with a leading gas turbine manufacturer. 

Aircraft gas turbines are subjected to rigorous tests 
immediately after manufacture, in order to ensure re¬ 
liable operation and give confidence in design predic- 
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tions. One such test is a dynamic verification test, con¬ 
ducted on a static test-bed such as that shown in Fig. 1. 
The dynamic models obtained can then be used for 
simulation and as a basis for control system design. 
Engine dynamics arise from complex, interacting phe¬ 
nomena: gas-flow behaviour in the compressor and tur¬ 
bine (affected by air inlet as well as engine conditions), 
shaft inertias and losses, fuel-flow transport delay, fuel 
dispersal and combustion, and the thermal behaviour of 
the engine and its surroundings. 

Identification of aircraft gas turbine dynamics has his¬ 
torically relied on the use of “wobble” tests, in which the 
engine is excited by single sines of different frequencies. 
The gain and phase shift are then calculated at each 
frequency, generating Bode plots which are easy to inter¬ 
pret and familiar to control engineers. Low-amplitude 
sinewave tests can be made insensitive to the nonlineari¬ 
ties in the engine response and to noise. However, they 
have severe drawbacks, the primary one being the ex¬ 
pense associated with very long test durations to allow 
the decay of initial transients at each frequency. 

Moreover, the presence of nonlinearities necessitates 
taking a frequency response at each of several operating 
conditions spanning the range of engine speeds from 
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Fig. 1. Aero gas turbine test-bed. 


ground idle to full power. There are at least two sources 
of nonlinearity in the engine response. The first of 
these is the nonlinear variation of the engine dynamics 
with operating power and condition (Saravanamuttoo, 
1992). Secondly, the dynamics change with the thermal 
state of the engine; for instance they differ between a fast 
and slow acceleration, due to effects such as changes in 
blade-tip clearances (Pilidis & Maccallum, 1986). Some 
such effects can be treated as slow linear modes but the 
frequency responses give very limited insight into the 
slow thermal dynamics of the engine. 

Finally, the current wobble-test procedure, although 
believed to give accurate estimation of the small-signal, 
steady-state engine behaviour, provides little indication 
of the accuracy of the results. 

The work described here is intended to improve the 
speed and quality of engine testing and to allow better 
exploitation of test-bed results. The overall aim is to gain 
insight into alternatives to traditional testing and identi¬ 
fication methods for multi-shaft aircraft gas turbine en¬ 
gines. This case study utilised a large set of engine records 
from a twin-shaft Rolls Royce Spey engine, logged by the 
Defence Evaluation and Research Agency (DERA) Pyes- 
tock, according to test schedules specified by the univer¬ 
sity teams. 

This paper concentrates on the dynamic relationship 
between the measured input fuel flow and the high- 
pressure (HP) and low-pressure (LP) shaft speeds, de¬ 
noted by N h and N L . The engine speed control was 
operated in open-loop and a perturbed fuel-demand sig¬ 
nal was fed to the fuel-feed system, which regulates the 
fuel flow to the engine by means of a stepper valve, as 
shown in Fig. 2. The fuel flow was measured downstream 
of the fuel-feed system, using a turbine flow meter, in 
order to exclude the fuel-feed dynamics from the engine 
model. 

Any alternative gas turbine identification technique 
must allow the user to assess the accuracy of the resulting 
models and must not yield lower accuracy than the 
current procedure. The intention is to improve on the 


“wobble” test technique currently used at Rolls Royce 
through: 

• Reduction in the time (and hence cost) of engine 
testing. 

• Improved accuracy of the estimated engine models 
and an assessment of the model uncertainty. 

• Greater insight into the slow engine dynamics. 

• Better understanding of the nonlinear behaviour. 

The combination of shorter, and thus cheaper, test 
runs and better coverage of engine behaviour is a strong 
incentive to look for more efficient test signal designs and 
model-estimation procedures. The three system identi¬ 
fication approaches addressed in this paper are: 

• multisine test signals and frequency-domain identifica¬ 
tion techniques, for improved linear modelling; 

• extended least squares with optimal-smoothing, for 
finding time-varying linear models; 

• multiobjective genetic programming, for the selection 
and identification of a nonlinear model structure. 

After a description of each approach and associated 
results, conclusions will be drawn about the scope and 
merits of the various techniques. 


2. Multisine signals and frequency-domain identification 

Previous work by Evans, Rees and Jones (1995) and 
Evans, Rees and Hill (1998) illustrated how multisine test 
signals and frequency-domain techniques could be used 
to accurately estimate parametric and nonparametric 
models of a gas turbine engine. That work focused on the 
study of the engine dynamics at a single operating point. 
In this section, the same techniques are used to estimate 
5-domain models of the engine dynamics at a range of 
points, with the aim of verifying the linearised thermo¬ 
dynamic models of the engine. 

These thermodynamic models are of great importance 
in the development stages of a gas turbine. Once the 
turbine is built, it is necessary to check the accuracy of 
the thermodynamic models with data from the actual 
engine. The models estimated using frequency-domain 
techniques have proved to be of great value in this task. 

2.1. Multisine testing 

It was important to reduce the time spent testing the 
engine and, to this end, multisine signals were used to 
excite the gas turbine. These are simply an arbitrary sum 
of harmonically related cosines 

F 

u(t)= X a k cos(i k 2nf 0 t + (j) k ) (1) 

k= 1 
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Fig. 2. Engine test facility. 


with a k the cosine amplitude, i k its harmonic number, 
f 0 the signal fundamental frequency, 4> k the harmonic 
phase and F the number of frequencies in the signal. The 
advantages of multisines over other test signals are 
numerous. The main benefits arise from the fact that 
multisines can be strictly bandlimited, which is very con¬ 
venient when no anti-alias filters can be used, and that 
the vector of harmonics can be arbitrarily defined, which 
allows the signal power to be placed exactly at the de¬ 
sired frequencies. In addition, a benefit common to all 
periodic signals is that the nonexcited spectral lines can 
be eliminated from the data set used for estimation. 

By including only odd harmonics in the signal (i.e. the 
fundamental, the third harmonic, fifth harmonic, and so 
on) the multisine is made immune to the influence of 
even-order nonlinearities. All the output frequencies gen¬ 
erated by even-order nonlinear terms, such as a squaring 
function, will fall at even harmonics and not affect the test 
frequencies (Evans et al., 1995). Such signals are termed 


odd-harmonic multisines. The presence of a weak even- 
order engine nonlinearity had been detected in a previous 
study (Evans et al., 1995). 

The harmonic phases must be carefully selected to 
minimise the crest factor (CF) of the signal 

max{|a(f)|} 

r.m.s.{u(0}’ 1 } 

which is a measure of the time-domain amplitude com¬ 
pression of the signal for a given root mean square (r.m.s.) 
value. Test amplitudes are often limited by safety or 
operational considerations (such as preventing surge or 
engine damage) and also by the need to avoid overly 
exciting the system nonlinearities. This translates to a set 
value of max{|w(f)|}. With such a limitation, a signal with 
a low CF will allow more power to be injected into the 
system and thus improve the signal-to-noise ratio (SNR). 
The CF minimisation technique proposed by Guillaume, 
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Schoukens, Pintelon and Kollar (1991) was employed, 
since it results in the lowest CFs achieved to date. 

The use of broad-band multisines allows a large num¬ 
ber of frequency points to be measured in one test, which 
is much faster than measuring frequency by frequency 
using single sines. This can be illustrated by considering 
the relative test time with each approach. For single 
sines, the total test time for a set of F equally spaced 
measurements, at integer multiples of a fundamental fre¬ 
quency /o, will be 


It was found that measuring six periods of the multisine 
was sufficient to achieve an acceptable accuracy, giving 
a test time of 

Multisine T test = ^6 + (8x 2.2) = 617.6 s. (7) 

This gives a relative test time of 

r 1 i r 

Relative T test = —— = 0.17, (8) 


T t est — M ss —F + T set F — (m ss — + T set ]F, (3) 

Jo \ Jo / 

where M ss is the number of averages and T set is the 
settling time allowed after each change of frequency 
before commencing measurements. This expression 
assumes that each sine is measured for the same 
time period, in order to ensure equal accuracy at each 
frequency. 

The total test time for a multisine signal will be 


1 

T — \/t _ T 

1 test ly± ms r ' 1 seti 

Jo 


( 4 ) 


where M ms is the number of averages of the multisine. 
Their relative test time can be expressed as 


Relative T test 


Multisine T test 
Single sine T test 


M ms T o + T set 
( M SS T 0 + T set )F 


where T 0 is the period of the fundamental frequency. 
It is assumed that averaging is performed over the basic 
interval T 0 in each test. It can be seen that the single sine 
tests will take F times longer than those with the multi¬ 
sines and that the relative benefit will increase with the 
number of test frequencies. However, this comparison 
takes no account of the relative accuracy of the tests. 
Since the power per harmonic is less for the multisine 
signal, it may be necessary to perform a greater number 
of averages in order to achieve an acceptable accuracy 
(M ms > M ss ). 

Illustrative example: Consider the design of a test to 
measure the dynamics of the LP shaft, having 5-domain 
poles at Si — — 0.45 and s 2 = — 1.8, with correspond¬ 
ing time constants t x = 2.2 s and t 2 = 0.55 s. In order 
to measure the FRF of this system a frequency range 
must be selected that adequately covers the system 
break-points. Since the break-points are at 0.07 and 
0.29 Hz then a frequency range of 0.01-0.6 Hz is seen to 
be adequate. A 30 odd-harmonic multisine with a funda¬ 
mental frequency of 0.01 Hz has a bandwidth of 
0.01-0.59 Hz. 

Selecting a settling time of eight times the slowest time 
constant, the minimum test time for single sines will be 

Single sine T test = + 8 x 2.2^ x 30 = 3528 s. (6) 


which shows that the use of a multisine signal will reduce 
the test time by greater than 80%. 

A range of multisine signals was used to excite the 
engine at different operating points. This included the 30 
odd-harmonic signal previously mentioned, with power 
between 0.1 and 0.6 Hz. Signals with power concentrated 
at low frequencies were also employed, in order to inves¬ 
tigate the presence of slow dynamics. The time records of 
one such signal are shown in Fig. 3, where four periods 
of the measured fuel flow and the HP and LP shaft 
speeds are plotted. The spectra of the measured fuel flow 
and the HP shaft speed are shown in Fig. 4. The input 
signal consists of odd harmonics between the funda¬ 
mental and the 71st harmonic, with increased spacing 
between the higher harmonics. The results are for a test 
at 75% of the maximum HP shaft speed (%N H ) with an 
input amplitude of + 10% of the steady-state fuel flow 
(%W f ). It can be seen that the SNRs are very good for 
this input amplitude. By excluding the noise lines and 
averaging over four periods, SNRs of greater than 40 dB 
were obtained for the input and output. 

By omitting all the even harmonics from the input 
signal, and also excluding some odd harmonics, it is 
possible to detect the influence of even- and odd-order 
nonlinearities by the presence of power at the excluded 
frequencies in the output (Evans, 1998). Any significant 
nonlinearities in the system will generate contributions at 
these excluded harmonics, which will rise above the noise 
floor. No significant components are visible in the output 
plot, which suggests that the nonlinear influence is small 
for this input amplitude range. 

Tests were conducted at a range of input amplitudes, 
from +1% W f to +20% W f . It was found that an 
amplitude of ± 10% represented a good compromise 
between the need to increase the amplitude in order to 
improve the SNRs and the need to limit the amplitude in 
order not to excite the engine nonlinearities. Such 
a trade-off must be made in testing most real systems. 

2.2. Frequency-domain identification 

The data obtained from testing with multisine signals 
were then used to estimate both nonparametric and para¬ 
metric models in the frequency domain. Since periodic 
signals were used, the frequency response function (FRF) 
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(a) 
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Time (s) 

Fig. 3. Time record of the measured (a) fuel flow, (b) FIP shaft speed and (c) LP shaft speed. 


(a) 



Frequency (Flz) 

Fig. 4. Spectra of the (a) measured fuel flow and (b) FIP shaft speed. 


was estimated as the ratio of the mean values of the 
output and input Fourier coefficients at the discrete test 
frequencies ca K 


H E v(jc>k) — 


(l/M)g.,r.(M) 

(l/M)S?„ CUM)’ 


where U m (ja > K ) and Y m (jco K ) are the input and output 
spectra measured across M periods of the input and 


output signals. This was termed the error-in-variables 
estimator by Guillaume (1992), who showed that it is 
both asymptotically unbiased and efficient in the pres¬ 
ence of normally distributed noise on the real and imagi¬ 
nary parts of the input and output Fourier coefficients. It 
has been shown by Schoukens, Guillaume and Pintelon 
(1993) that the variance of the estimated FRF is inversely 
proportional to the number of measurements and the 
power of the input harmonics, and proportional to the 
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(a) 



Freq (Hz) 
(b) 



Freq (Hz) 


Fig. 5. FRFs for the HP shaft at 65% N H (o), 75% N H (x) and 85% N H ( + ). 


(a) 



Freq (Hz) 
(b) 



Freq (Hz) 


Fig. 6. FRFs for the LP shaft at 65% N H (o), 75% N H (x) and 85% N H ( + ). 


noise variances referred to the system output. With such 
high SNRs, of greater than 40 dB after averaging, the 
uncertainty of the FRF estimates was very small. 

Fig. 5 shows the HP shaft FRFs resulting from tests at 
three different operating points, which illustrates the 
evolution of the shaft dynamics with the shaft speed. 
Fig. 6 shows the equivalent evolution for the LP shaft. 
The FRFs give a good insight into the engine dynamics 


and they show that both shafts have steady-state gains 
that decrease, and dynamics that become faster, as the 
shaft speeds increase. 

Parametric frequency-domain identification involves 
selecting the parameters of an s-domain model: 


H(s) 


bo + b^s + ••• + b m s m ^- sTd 
ciq T ci^s T ••• T a n s 


( 10 ) 
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Fig. 7. Evolution of the estimated poles (x) and zeros (o) with operating point: (a) HP shaft and (b) LP shaft. 



75 

Shaft Speed (%NH) 



Shaft Speed (%NH) 


where m is the number of zeros, n the number of poles 
and T d the pure time delay. The estimator developed by 
Pintelon, Guillaume, Rolain and Verbeyst (1992) and 
implemented by Kollar (1997) was used for this purpose. 
The estimator has a least-squares form, with each fre¬ 
quency point weighted by a noise variance term. The 
pure time delay T d can also be included as a free para¬ 
meter for estimation, which is an attractive feature of the 
frequency-domain approach, since its value is not fixed to 
multiples of the sampling interval. 

A nonparametric noise model is employed and the 
noise variances and covariance are required as a priori 
information. These can be estimated from the data before 
estimation, provided several independent measurements 
are available. It has been shown that the estimated mod¬ 
els are strongly consistent, for an increasing number of 
data points in each experiment, if the number of indepen¬ 
dent measurements is greater than or equal to four 
(Schoukens, Pintelon, Vandersteen and Guillaume, 
1997). In the gas turbine testing, each period of the test 
signal was considered as an independent measurement 
and at least four periods were measured in each test. 

Tests were conducted at several operating points along 
the turbine running range from 53 to 90% N H and para¬ 
metric models of the HP and LP shafts estimated at each 
point. All of the estimated models had real and negative 
poles and zeros. The model selection and validation 
procedure is described in detail in Evans, Rees and Bor- 
rell (2000). The locations of the poles and zeros are 
presented in Fig. 7 and some clear features can be de¬ 
duced from the plot. 

Consider first the tests at 53 and 55% N H , where the 
locations of the poles and the zero suggest that the 


dynamics are similar for both the HP and LP shafts. At 
the operating points 65 and 75% N H the plots suggest 
a first-order model for the HP shaft dynamics and a sec¬ 
ond-order model for the LP shaft dynamics. An addi¬ 
tional pole-zero pair was found at low frequencies in 
each model. These low-frequency pole-zero pairs were 
also found by Evans (1998) when testing the turbine at 
75% N h , where it was suggested that they were modelling 
a slow dynamic such as heat soakage. Indeed, the time 
constants associated with this pair are between 5 and 
22 s, which are too slow to be modelling the shaft dynam¬ 
ics at these operating points. 

At the higher operating points, 85-90% N H , a different 
trend is observed for the HP shaft. The low frequency 
pole and zero have moved apart making the two poles 
more distinct. In fact, they suggest that the HP shaft 
behaves like a second-order system at these operating 
points, as it does at low speeds. This may express the fact 
that the interaction between the shafts is greater in these 
regions. The LP shaft is also second order at 85 and 90% 
N h but with significantly different dynamics to that of the 
HP shaft. 

2.3. Verification of thermodynamic models 

Engine models are required both in the development 
and operational stages of the life of a gas turbine. Ther¬ 
modynamic models are derived during the development 
stage, based on knowledge of the engine physics, and 
provide important insights into the engine behaviour. 
Such models are both complex and nonlinear, making 
them unsuitable for use in the design of engine control 
systems. 
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Table 1 

Modes of full-order linearised thermodynamic model 


S-plane value 

hi 

Frequency 

(Hz) 

Time constant 
(ms) 

0.69 

0.11 

1451 

3.1 

0.50 

319 

61 

9.7 

16 

107 

17 

9 

335 

53 

3 

403 

64 

2 


It is therefore common practice to linearise the thermo¬ 
dynamic models around a series of operating points and 
then carry out a model-order reduction, in order to arrive 
at models which are suitable for control system design. 
Since these models are based on a priori assumptions 
about the engine physics it is important to then validate 
their performance against real engine data. 

In work conducted for Rolls Royce pic, Jackson (1988) 
showed that, for a given stationary operating point, the 
higher-order nonlinear thermodynamic models can be 
reduced to linear models with the same order as the 
number of engine shafts. The models are first linearised 
using small perturbations, which in the case of the Spey 
engine results in models with five inputs, 10 outputs and 
15 states. Table 1 shows the four slowest and two fastest 
modes of such a full-order model, linearised at an operat¬ 
ing point of 75% N h . 


There is a big difference between the time constants 
of the two slowest modes, associated with the shaft 
dynamics, and the other modes, which are associated 
with thermodynamic processes in the gas stream. 
A model reduction procedure can thus be employed to 
eliminate all of the faster modes. A complete library of 
such models was generated for the Spey engine, across 
a range of operating points, by staff at section APD5 
(1993) of DERA Pyestock. These linearised, reduced- 
order models will be referred to simply as the thermodyn¬ 
amic models for the remainder of this paper. Evaluating 
the transfer function matrices of the state space models 
allows the HP and LP shaft dynamics to be expressed in 
transfer function form, with common poles but different 
zeros. 

A comparison can now be made between the thermo¬ 
dynamic models and the models estimated directly from 
engine data. The poles and zeros of the thermodynamic 
models and the estimated models are plotted together 
in Fig. 8. In the case of the HP shaft thermodynamic 
models, the second pole and the zero lie over each other 
in the range 55-75% N H . The low-frequency pole-zero 
pairs of the estimated models have been omitted from 
this comparison, since they are not related to the shaft 
dynamics. The plot shows a number of interesting trends. 
In the first place, there is good agreement between the 
thermodynamic models and the estimated models with 
regard to the first-order dynamics at low shaft speeds. 
This is true for the pole positions of both the HP and 
LP shafts. A deviation begins to appear at the higher 
operating points. 


(a) 




Fig. 8. Poles and zeros of thermodynamic and estimated models for the: (a) FIP and (b) LP shafts. Thermodynamic models: poles (solid) and zeros 
(dashed), and estimated models: poles (x) and zeros (o). 
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The most notable discrepancy is found in the loca¬ 
tion of the second pole and the zero of the LP shaft. 
The thermodynamic models are clearly not model¬ 
ling the evolution of this pole and zero with shaft 
speed in the way suggested by the estimated models. The 
estimated models are based on real data gathered 
from an engine and the models show close agreement 
with the engine FRFs estimated at each operating 
point. It can thus be argued that the estimated models 
reflect the true dynamics of the engine and that the 
thermodynamic models, which are based on a priori 
assumptions about the engine, are problematic. Greater 
confidence should thus be placed in the experi¬ 
mental results, which clearly show that the thermo¬ 
dynamic models are not adequately representing the 
dynamics of either shaft at higher operating points. 
A particular disagreement has been found for the sec¬ 
ond-order dynamics of the LP shaft, which are of special 
interest since they model the interaction between the 
shafts. 

It is also interesting to note the behaviour of the 
estimated HP shaft models at low and high operating 
points. They suggest that the HP shaft exhibits 
second-order dynamics in these regions. This second- 
order behaviour is predicted by the thermodynamic 
models at the higher operating points but not at the 
lower. 


3. Extended least-squares with optimal smoothing 

Traditional single sine engine testing deals with non¬ 
linearity only by generating small-signal models about 
a range of steady-state operating points. It does not cover 
effects far from equilibrium, which may be prominent in 
large transients such as slam accelerations (rapid, full 
throttle opening). Small-signal models for short-term 
response also give little insight into the slow thermal 
dynamics of the engine. Although such models are valu¬ 
able for verifying engine performance near steady state, 
they are a poor basis for control-system design for large 
inputs. The main aim here is to investigate large-transient 
behaviour and slow thermal dynamics and to see whether 
they can be identified in black-box models from relatively 
simple large-perturbation tests, economical in test 
duration. 

All records were digitally low-pass filtered with 
a bandwidth of 1.6 Hz, removing substantial wide¬ 
band noise from the flow records, and resampled at an 
interval of T — 0.1 s. The sample means were removed to 
avoid the need for an offset term. Initial fitting by ex¬ 
tended least squares of first- and second-order time- 
invariant discrete transfer-function models to results of 
inverse-repeat maximum-length binary sequence 
(IRMLBS) tests with small ( ± 10% W f ) perturbation 


about high-pressure shaft speeds from 53 to 90% 
N h showed: 

• poor observability of the small, HP shaft second 
mode, but a larger, better identified second mode 
for the LP shaft. Biases are in the order of 10 and 
25% in HP and LP dominant time constant if the 
second mode is omitted, even though that mode is 
around 10 times faster and has steady-state gain typi¬ 
cally 30 times smaller for the HP shaft (although only 
3-4 times smaller for the LP shaft at the highest 
speeds) 

• over most of the speed range, lower small-signal 
dc gain from the transfer-function models than from 
linearizing the static fuel-flow vs. speed characteristic. 
Differences are around 10 and 10-20% for HP and 
LP shafts respectively. The reverse applies at each end 
of the speed range, by 20-30%. This all suggests that 
the MLBS tests fail to elicit all the low-frequency 
response 

• smooth variation of the parameters with speed, and 
low standard errors for the parameters of the domi¬ 
nant mode in second-order models. 

For large-transient tests, an optimal-smoothing gener¬ 
alisation of extended least squares (Norton, 1975, 1986) 
yields time-varying transfer-function models of approx¬ 
imately uniform quality throughout each test. The idea is 
that the time variation reveals nonlinearity and slow (e.g. 
thermal) dynamics that would be averaged out by 
a time-invariant, linear, low-order model. Space only 
allows reporting of results for N H . 

First, a 254-bit IRMLBS varying N H by + 10% was 
superimposed on a triangular wave of the same period, 
100 s. This drives N H between 65 and 85% after several 
minutes’ initial “soaking” at 65% N H . The most striking 
feature of the estimates in Figs. 9 and 10 for the pole 
p and steady-state gain g of the dominant mode is the 
initial rise in p, rapid during the initial soak (of which 60 s 



Fig. 9. Dominant discrete pole (p), IRMLBS + triangular wave. 
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Fig. 10. Dominant mode gain ( g ), IRMLBS + triangular wave. 



Fig. 12. Dominant mode gain ( g ), triangular wave only. 



Fig. 11. Dominant discrete pole ( p ), triangular wave only. 


is shown) then slower. It is explicable as a slow res¬ 
ponse to cooling of the engine during the 5-min 
soak following a test at speeds between 75 and 91%, and 
implies a long thermal time constant. By contrast, 
g is near-constant at about the value expected for 
65% N h while soaking, then varies about the value ex¬ 
pected at 75% N h , the mean speed. However, the peri¬ 
odic variation is less than expected from the small-signal 
results, indicating that thermal equilibrium is never 
reached. 

Another unexpected feature is that the variation of 
g and, less markedly, p includes higher-frequency compo¬ 
nents in each perturbation cycle. A possible reason is 
local structure within the IRMLBS sequence, so com¬ 
parison with results with no IRMLBS imposed on the 
triangular wave, Figs. 11 and 12, is called for. Figs. 11 and 
12 show smoother variation than Figs. 9 and 10 but 
similar modulation effects. Clearly dynamics which can¬ 
not be represented by moving through a succession of 
low-order, small-signal models are present. They may be 
ascribed to higher-order dynamics, perhaps requiring 
second and/or higher derivatives of speed in the differen¬ 


tial equation describing shaft-speed response to fuel flow, 
e.g. in the torque equation 


power m 


N, 


= JN h + B(N h )N h 


H 


+ 


power to exhaust + heat loss 


N 


H 


(ID 


Such dynamics are ignored in the small-signal models. 

A final conclusion to be drawn from Figs. 9 and 11 is 
that dynamics with a time constant of about 150-200 s 
are present, evidenced strongly by the overall trend of the 
pole. These dynamics have a fairly small effect on g but 
influence the dominant fast-response time constant 
— T/\np by about 30% in Fig. 11. 


4. Nonlinear system identification using multiobjective 
genetic programming 

Genetic programming (GP) is an evolutionary para¬ 
digm where the computer structures which undergo 
adaptation are themselves represented as computer pro¬ 
grams (Koza, 1992). Koza immediately recognised that 
symbolic regression was as an obvious application do¬ 
main, where mathematical expressions and their numeric 
coefficients may be obtained to provide a good fit to a set 
of data points. Rodriguez-Vazquez and Fleming (1997) 
describe how GP can be used to search for a suitable 
nonlinear model structure for a NARMAX model 
(Leontaritis and Billings, 1985) of a dynamic system, 
given an input and output data set. Here, the multiobjec¬ 
tive genetic programming (MOGP) approach of 
Rodriguez-Vazquez and Fleming (1998) is used to 
obtain NARMAX models to attempt to describe the 
global behaviour of the gas turbine engine over its 
operating range. The MOGP fitness function, defined 
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Table 2 

MOGP system identification objectives 


Attribute 

Objective 

Description 

Model complexity 

1. Model size, NT 

2. Degree of nonlinearity, DG 

3. Model lag, LG 

Number of process and noise terms 

Maximum order term 

Maximum lagged input, output and noise terms 

Model performance 

4. Residual variance, VAR,- 

5. Long-term prediction error, LTPE ; 

Variance of the predictive error between the predicted and the 
measured outputs for dataset i 

Variance of the long-term prediction error for dataset i 

Model validation 

6. ACF 

7. CCF 

8-10 HOCj as defined in Eq. (16) 

Auto-correlation, cross-correlation and higher-order correla¬ 
tion-based functions 


as a multiobjective function, generates a set of non- 
dominated models and, thereby, a set of alternative 
system models. The complexity, performance and also 
statistical validation of the models are computed for each 
MOGP population member and contribute to their fit¬ 
ness evaluation. The final set of nondominated model 
solutions is able to capture the dynamic characteristics 
of the system and reproduce its behaviour at different 
operating conditions. 

The model representation used is the NARMAX 
model of Leontaritis and Billings (1985). The general 
NARMAX model is defined as a nonlinear function of 
the input, output and noise signal terms as given by the 
following equation: 

y(k) = F <f (y(k - 1), ...,y(k — n y ), u{k - 1), ..., u{k - n u ), 

e{k - 1), ...,e{k- n c )) + e{k), (12) 

where y(k), u(k ) and e(k) represent the output, input and 
noise signals, respectively, n y , n u , and n e are their asso¬ 
ciated maximum lags, and F is some unknown nonlinear 
function of degree /. 

Based upon the model representation defined by 
Eq. (12), the MOGP population consists of tree- 
structured individuals that readily represent alternative 
structures for the application of the NARMAX ap¬ 
proach. Potential models are encoded as hierarchical 
tree structures, thus providing a dynamic and variable 
representation, and these constitute members of a popu¬ 
lation of different model structures. These structures con¬ 
sist of functions (internal nodes) and terminals (leaf 
nodes) appropriate to the problem domain. Hence, the 
function set is here defined as F — {ADD, MULT} = 
{+,*}, and the terminal set as T = (X 0 , ...X ny , 
X ny +15 X ny + nu X ny + nu -|- 1 ,..., X ny + nu 4 . ne } \ c,y(k 1), 
...,y(k- n y ), u(k - 1), ...,u(k - n u ), e(k - 1), ...,e(k- n e )}. 
For example, one hierarchical tree representation 


of the polynomial NARMAX model might be expressed 
in Polish notation as (ADD(ADD(ADD XIX4) 
X5)(MULT (ADD X2 X3)(ADDA 1X2))). This is equiva¬ 
lent to the polynomial nonlinear model defined as 

y(k) — 6 0 + 6 1 y(k — 1) + d 2 y(k — 2) + d 3 u(k — 1) 

+ 0 4 e{k - 1) + 6 5 y(k - l) 2 

+ 0 6 y(k - 1 )y(k - 2) (13) 

where {Xl,X2,X3,X4,X5j = (1.0(the constant term), 
y(k — 1 ),y(k — 2), u(k — 1 ),e(k — 1)}. A least-squares 
algorithm is applied to compute the parameters, 0j ; 
to minimise the residual of errors, s, between the 
measured output y(k) and the predicted output y(k) that 
is given by 

s(k) = y(k) - y(k,§). (14) 

In order to perform selection at each generation in the 
MOGP approach, all model measures considered in the 
identification are evaluated for each member of the popu¬ 
lation. The fitness value of each population member is 
assigned by means of a rank-based fitness method (Fon¬ 
seca and Fleming, 1998). This fitness evaluation is based 
on the definition of Pareto-optimality or nondominance. 
If we consider a minimisation problem and, given two 
71 -component objective function vectors, f„ and f v , we can 
say that f„ dominates f, (is Pareto-optimal) if 

Vie(1, f Ui <f Vi a 3ie{l, ...n}, f Ui <f Vi , (15) 

producing a set of possible and valid solutions known as 
the Pareto-optimal or nondominated set. That is, a Par¬ 
eto-optimal solution is one in which an improvement in 
one of the design objectives will lead to degradation in 
one or more of the remaining objectives. 

Thus, the model complexity, performance and vali¬ 
dation attributes of the modelling process can 
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simultaneously be evaluated by considering them in the 
multiobjective fitness function. The objectives used in 
this study are defined and classified as shown in Table 2. 
The validation attributes are measured, based upon the 
auto-correlation (ACF), cross-correlation (CCF) and 
higher-order correlation (HOC) functions proposed by 
(Billings & Voon, 1986) and defined as 

ACF = <P ee (x) = £[e(k — t)s(/c)] = <5(t), 

CCF = <2> M£ (t) = Elu(k - x)e (k)] = 0 Vt, 

HOCi = (P £(£U) (t) = E[&(k)&(k — x)u(k — t)] = 0, t > 0, 

HOC 2 = <P l g E {x) = E[{u\k ) - £[u 2 (k)]}fi(k - t)] = 0 Vt, 

HOC 3 = <2>„ v (t) = E[_{u 2 {k) - E[u 2 {k)~\}&\k - t)] = 0 Vt. 

(16) 

The objectives related to the performance and 
complexity of candidate NARMAX models are shown in 
Table 2. 

The identification approach presented here concerns 
the identification of nonlinear models that can cope with 
a specific range of system operation. As such, it is based 
upon data obtained from multisine excitation signals at 
different operating conditions. The objective is to find 
a set of NARMAX models that minimises the n-compo- 
nents vector criterion which includes measures of vari¬ 
ance and LTPE at these different operating conditions, 
albeit, operating within a specific working region. 
The multiobjective vector function to be minimised is 
defined as: 

f = [NT, DG, LG, VAR,, LTPE,, ACF, CCF, HOC 7 ] T 

where the variables are defined in Table 2. 

Initial runs that were conducted with a multiobjective 
vector function 

f = [NT, DG, LG, VAR,, LTPE,] t 

which did not include model validation criteria, produc¬ 
ed a set of linear models. However, these solutions 
failed to satisfy validation requirements. Consideration 
of the full objective function vector produced models 
containing quadratic terms. Since multiobjective opti¬ 
misation produces a family of solutions, the MOGP 
procedure invites designer interaction. This resulted in 
three candidate solutions, whose term content is given 
in Table 3. Note that even though the terminal set T 
included past values of the residuals {e(k — 1), ..., 
e(k — n e )j, no NARMAX structures emerged; instead, 
only NARX model structures were produced. These three 
models were subjected to rigorous simulation tests 
(Rodriguez-Vazquez 1999) and Model 3 was selected as 
the most appropriate. The structure and parameters of 


Table 3 

Quadratic polynomial models 


Model term 

Model 



1 

2 

3 

c 

▲ 

A 

A 

y(k - 1) 

▲ 

A 

A 

y(k - 2) 

A 

A 

A 

u(k — 1) 



A 

u(k — 2) 



A 

y(k- l) 2 



A 

y(k - 2) 2 

A 



y(k — 1 )u(k — 1) 

A 

A 


y(k — 2 )u(k — 2) 

A 

A 



Three-level periodic sequence, 58-65-70% NH 

751-i-i-i-i- 



55 1 - ‘ -i-1-1-1-1 

0 20 40 60 80 100 120 

Residuals 

1 i 



-0.5 1 - 1 - 1 - 1 - 1 - 

0 20 40 60 80 100 120 

Time [s] 

Fig. 13. (a) Three-level periodic sequence of period 100 s, each period 
consisting of quarter-periods at levels 65, 70, 65 and 58% N H (actual: 
solid line, estimated: dotted line), (b) Corresponding residuals. 


this quadratic NARX engine model are 

y(k) = - 1.16236 x 10° + 7.96090 x 10“ - 1) 

+ 2.39363 x 10“ V(k-2) 

+ 4.33425 x 10“ 3 u(k - 1) 

- 5.90401 x 10“ 4 u(k - 2) 

x 4.37337 xl0“ 4 y(k-l) 2 . (17) 

The ability of this model to represent the engine relation¬ 
ship between the measured fuel flow and the high pres¬ 
sure shaft speed on two exacting test scenarios is 
illustrated in Fig. 13 (three-level periodic signals over 
the range 58-65-70% N H ) and Fig. 14 (a large transient 
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Fast accel-steady-fast decel, 55-85% NH 



Fig. 14. Fast acceleration/deceleration signal (55-85% N H ), actual: 
solid line, estimated: dotted line. 


55-85% N h ). The upper plot in Fig. 13 shows the 
two plots of actual and estimated data (nearly sup¬ 
erimposed) and the lower plot displays the residuals 
(difference between the actual and estimated data) 
which are very small. Recall that the objectives directed 
the search to find good low-order simple non-linear mod¬ 
els. Fuller results are provided in Rodriguez-Vazquez 
(1999). 


5. Conclusions 

The three techniques outlined in this paper were all 
applied to real engine test data. They provide significant 
insights into alternative identification strategies. The 
aims, to improve the efficiency of gas turbine dynamic 
testing and to see how far the identification approaches 
can give insight into complex engine behaviour, have 
largely been realised. The conclusions from each of the 
techniques will be summarised in turn. 

The contribution based on the use of multisine signals 
and frequency-domain techniques has been: 

• The design of multisine signals with a low crest factor 
to excite the dynamics of the gas turbine. These multi¬ 
sines allow a very significant reduction in test time, 
compared with the single-sine “wobble” tests currently 
used, without any significant loss of accuracy. Reduc¬ 
tions in test times of 80%, or greater, can be achieved 
by using multisine signals. 

• The levels of noise and nonlinearities present in the 
turbine have been assessed for small-signal tests. For 
an input signal amplitude of + 10% W f , it was found 
that the SNRs were better than 40 dB, after averaging, 
and that the nonlinear effects were negligible. 

• The nonparametric frequency response functions have 
been estimated at various operating points and their 
confidence bounds have been determined. These were 
found to be very small for an input signal amplitude of 


± 10% Wf, providing a high degree of confidence in 
the estimated frequency responses. Parametric 5-do- 
main models have been estimated directly in the fre¬ 
quency domain and the uncertainty on their poles and 
zeros was found to be small. The pure time delay was 
included as an estimated parameter and the results 
matched those predicted by theory (see Evans et al., 
2000 ). 

• The parametric models estimated from the engine data 
were then used to verify the linearised thermodynamic 
models derived from the engine physics. This showed 
significant discrepancies between the estimated and 
derived models, which warrant further investigation. 
This was particularly true for the second-order dy¬ 
namics of the LP shaft and the HP shaft at the lower 
and higher operating points. 

The extended least squares with optimal-smoothing 
approach allowed investigation of engine dynamics over 
large perturbations. Hitherto, gas turbine identification 
techniques have concentrated on capturing dynamics 
using small-signal models. Evidence of slow thermal 
dynamics, with implications for engine test duration, was 
found using this approach. 

In contrast to the second technique, which was con¬ 
cerned with time-varying linear models, the last tech¬ 
nique employed multiobjective genetic programming 
for the identification of time-invariant nonlinear 
(NARMAX) models. The multiobjective approach per¬ 
mitted the simultaneous evaluation of multiple indepen¬ 
dent criteria in the identification procedure and genetic 
programming was used to conduct the search over 
a range of nonlinear terms. The use of model complexity 
objectives alongside model performance objectives sup¬ 
ported the identification of low-order nonlinear models 
for a specific large-amplitude operating region. The vali¬ 
dation functions (correlation tests) led to the identifica¬ 
tion of second-order dynamic models containing 
a quadratic nonlinear term. By enabling designer interac¬ 
tion in the process, it was possible to selectively reduce 
the search space while pursuing a set of objectives. This 
focused the genetic search into a region of acceptable 
models. Validation of the identified models against differ¬ 
ing data sets proved the acceptability of the identified 
models. 
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